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Abstract. Observations of HST and groundbased data strongly suggest that most galaxies harbour central su- 
permassive black holes and that most galaxies merge with others. Consequently a black hole binary emerges as 
the two black holes are spiraling into the center towards each other. In our work we are investigating two ba- 
sic questions of our understanding of the central activity of galaxies and find that both can be answered with 
"yes": (1) Do the black holes actually merge? (2) And does the effect of the torque of the black hole binary on 
the surrounding stellar distribution help to explain the presence of the ubiquitous torus of molecular material 
surrounding apparently all active galactic nuclei? The first question is the topic of the present article, while the 
second question will be subject of a subsequent paper. Simulating the evolution of a stellar cluster in the potential 
of such a binary by solving the equations of the restricted three body problem we obtained the following results: 
Provided that the cluster is about as massive as the black hole binary the two black holes coalesce after ~ 10 7 yr 
due to ejection of stars and finally via emission of gravitational radiation. Whether a star is ejected or not crucially 
depends on its angular momentum. Almost all stars whose angular momentum is twice as large as that of a star 
circulating around the binary in a distance corresponding to that between the black holes, stay bound to the 
binary. In a sequence of models where the mass of the secondary black hole increases while Mi is fixed, a bigger 
fraction of stars is ejected. For a more massive binary also the cluster has to be more massive in order to allow 
the two black holes to coalesce. The merger then proceeds on smaller time scales. The cluster is depleted in the 
central region and the final distribution of stars assumes a torus-like structure, peaking at three times the initial 
distance of the two black holes. The relationship of the bound stars to the obscuring torus in active galactic nuclei 
will be investigated in a subsequent paper. 
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1. Introduction 

The bright radiation in the centers of galaxies is gener- 
ally thought to be produced by a powerful engine, which 
is loca pcd at the dynamical center of its host. Because 
the region from where this high luminosity is emitted is 
concentrated to a small volume, non-stellar activity is sus- 
pected to cause such a great energy output. It is argued 
that these active galactic nuclei (AGN) are powered by 
accretion onto a massive black hole (BH) in the center. 
While matter is sinking down in the potential of such a 
black hole, energy is released and this accretion process 
is thought to be dominated by an accretion disk. The 
maximum mass of these central objects, of order 10 -25 
times the stellar spheroidal mass of their host galaxies 



(Magorrian et al. 1998; 


Richstone et al. 1998 




Wang & 


Biermann 1998 




Macchetto 199E), is consistent with the 



mass in black holes needed to produce the observed en- 
ergy density in quasar light if reasonable assumptions are 
made about the efficiency of the energy production of the 



nucleus (Novikov & Thorne 1972; Shakura & Sunyaev 



1973; Chokshi fc Turner 1992; Blandford 1999). One of 



the best candidates of galaxies showing strong evidence 
for harbouring a massive BH in its center is NGC 4258. 
From VLBI observations of megamasers in its nucleus, 
Miyoshi et al. (1995Q deduce the existence of a mass of 



about 3.6 x lO^M^ which must be located inside the inner 



radius of ~ 13 pc. Peterson fc Wandcl (1999 ) use emission- 
line variability data on NGC 5548, a Seyfert 1 galaxy, to 
infer the existence of a mass of order 7 x 1O 7 M0 within 
the inner few light days of the nucleus. There are numer- 
ous other examples suggesting that most galaxies keep a 
super-massive black hole in their center. This is the first 
of our two basic assumptions. 
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According to hierarchical models a typical bright 
galaxy has merged a few times since the quasar era ( Frcnk 
et al 



98? 



Carlberg & Couchman 198E ; Baugh, Cole, & 



Frenk 996; Richstonc ct al. 1998). For binary black holes 
with M > 10 7 Mq the timescales for a decay in the dif- 



big blue bump (BBB) and the broad emission line region 
(BLR) as well as the narrow emission line region (NLR) 
can be seen directly On the other hand, objects classified 
as Type 2, show strongly absorbed X-rays and no clear 
evidence for the big blue bump and the BLR (Lawrence 



ferent regimes indicate that the binary will merge on a 



scale short compared to the time when the next merger 
event occurs (Xu & Ostriker 1994). The rate of mergers 
for bright galaxies in the observable universe may exceed 
1/yr ( Richstone et al. 1998] ) . Taking into account that lo- 
cal quasars might have been form ed by major mergers 
betwee n galaxies, Taniguchi (1999 ) proposes that all ac- 



& Elvis 1982). But the NLR is indistinguishable from that 
detected in Type 1 objects. This region is located at big- 
ger radii from the center and thus not blocked by matter 
at a smaller distance in the line of sight. Sometimes a 
Type 1 spectrum is seen in the polarized light of Type 2 
classified objects and thus gives strong evidence for the 



existence of an obscuring torus, see e.g. Miller, Goodrich 



tive ga |lactic nuclei of the local universe are products ei- fc Mathews (1991|) : this is interpreted within the unifica- 
tion model as radiation emerging from the center and then 
scattered by ionized matter in the opening of the torus 
into the line of sight (|Antonucci fc Miller 1985| ; |Miller fc 
Goodrich 19~90] |Tran fe Kay 1992j ). This model ascribes 



ther from minor or major mergers. This appears to be 
consistent with almost all important observ ational prop 
crtics of Seyfert galaxies. Recent studies (Mulchacy fc 



Regan 1997; Ho, Filippcnko, fc Sargent 1997; dc Robertis 



Yee, fc Hayhoe 1998 ) show that Seyfert galaxies do not 
have always either bar-structure or companion galaxies, 
what often had been considered to be the cause of effec- 



tive gas fuelling to the central engine (see, e.g., Adams 
1977], |f oguchi 1988|, jBhlosman, Begclman, fc Frank 1990j 



Barnes fc Hernquist 1992] ). An alternative way to supply 
the central engine with matter is the merger driven fu- 



tile different appearance of AGN rather to different orien- 
tations than to intrinsic differences. Its basic ingredients 
which are responsible for the non-spherical symmetry of 
the nucleus are relativistic beaming in a jet, the obscura- 
tion of the central engine by a molecular, dusty torus, and 
possibly the power of the central engine (see e.g the re- 
views by Antonucci 1993; Urry & Padovani 1995; Madcjski 



, see for example Toomre & Toomre (1972 


), 


Roos 


1998; 


Wills 1999 and articles by F 


alcke & Biermann 1995; 


) Mihos & Hernquist (1994), dc Robertis, Yee, & 


Falcke, Malkan, & Biermann 1995 


; Laniguchi 199G|; Rudgc 



(1981) 



Hayho ; (1998|). In order to complete a minor merger about 
10 9 yr are needed. This is enough time to smear out the 
relics so that most of the advanced minor mergers may 
be observed as ordinary- looking isolated galaxies (Walker, 
Mihos, & Hernquist 1996) and Seyfert galaxies are ob- 
served to possess a statistically signi ficant excess of faint 
(My > —1 8) companion galaxies ( Fuentes- Williams fc 



Stockc |l988[ ). The frequent merging of galaxies is our sec- 
ond basic assumption and consequently leads to massive 
black hole binaries in combination with the first assump- 
tion. 

Observations show that all ellipticals have central den- 
sity cusps which can be divided into 'weak' cusps (p cx 
r -(o. 3-i.i )^ bright ellipticals) and 'strong' (p oc r~ 2 , fainter 



radio-quiet quasar Mrk 205 might provide the first ev- 
idence for a dust-torus in a quasar (Reeves et al. 2001). 
XMM-data of IRAS 13349+2438 suggest column-densitie s 
which are consistent with a dusty torus (Sako et al. 2001). 
Observations impose important constraints on the prop- 
erties of such tori. From X-ray absorpt ion a column den- 
sity of N H w 10 24 cm~ 2 is deduced ([Mushotzky 1982; 



Mulchaey, Myshotzky, fc Weaver 1992j ; |Madejski 1998|) , 
i.e. the tori have to be optically thick. Number ratios of 
Type 1 and 2 ob jects tell us that such tori are also geomet- 
rically thick (see e.g. Taniguchi fc Anabuki (1999 ) and ref- 
erences therein) . The AGN spectra also show a prominent 



elliptic als) cusps ( Kormendy et al. 1994 Kormendy et al 



infrared bump between 1mm and 2pm (Sanders et al 



1996 ]_ [ 'aber et al. 1997 ). It has alread y been suggested by 
Begclman, Blandford, fc Rees (19~80| ) that the mass ejec- 
tion due the merging of two massive black holes should 
reduce the central density and that the core expands. In 
several numerical simulations it has been found that the 



merging of two galaxies and their central black holes can 
account for the weak cusps, where the sinking of the BH 
is the dominant mechanism (Makino & Ebisuzaki 1996; 



1989 ), which is thought to be generated by dust repro- 
cessing the incident radiation from the center. In about 
1 pc distance the dust in the torus is heated to its evap- 
oration temperature of about 1500 K and then reemitting 
the central ratiation in the infrared ( Banders et al. 1989|; 
Haas et al. 2000| ). |Pier fc Krolik (1992| ) and |Picr fc Kroht 



Nakano & Makino 1999; Milosavljevic & Merritt 2001 



Because the critical central regions of AGN are 
strongly nonspherical but spatially unresolved, orientation 
effects are thought to be important, sec e.g. the discus- 
sion of quasar spectra by Sanders et al. (198S ) . Thus the 
same object seen from other angles would be classified 
differently. In Type 1 objects the radiation emerging from 
the center is not blocked by matter in the line of sight 
of the observer. Hence the X-rays in the keV-range, the 



(1993) find that a very compact and thick dust torus, 
supported by radiation pressure, is not able to provide 
the required broad temperature range of the dust in or- 
der to explain the infrared (IR) emission. But a clumpy 
torus, as suggested by Krolik & Begelman 1988, has the 
advantage that the dust, when organized in clouds, can 
survive more easily the strong radiation field of the AGN. 
Also the dust in the outer parts can be heated to higher 
temperatu res and a broader temperature rang e can be 
achieved ( Efstathiou fc Rowan- Robinson 1995 ). Such a 
clumpy torus model seems to b e the most promising sce- 
nario also to Haas et al. (2000 ). A torus is also used in 
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high energy physics according to Protheroe & Biermann 
(1997^ who show for Blazars that the TcV 7-emission site 
is far from the central engine due to 77-interactions with 
the far infrared photon field of the torus. 

Now, based on our two underlying assumptions, nec- 
essarily leading to a black hole binary, we here investigate 
the influence of a stellar cluster and a black hole binary 
(BHB) in its center on each other in order to answer two 
basic questions in our understanding of the central activ- 
ity of galaxies: 

— Do the two black holes in the center of two colliding 
galaxies actually merge? 

— Does the effect of the torque of the binary black holes 
on the stellar distribution help to explain the presence 
of the ubiquituous torus of molecular material sur- 
rounding apparently all active galactic nuclei (AGN) 
in the "unified scheme"? 



Rces (1994) find the eccentricity never to exceed a value 
of order 0.1. After the stars surrounding the secondary 
black hole have been stripped off by tidal effects, it moves 
on bound orbits in the mean gravitational field of the 
core of the primary galaxy. Dynamical friction extracts ki- 



netic energy from the secondary BH ( Chandrasekhar 1943 



Binncy & Tremainc 1987) which binds to Mi at the radius 
7'b = (Mi/A/ corc ) i/a r corc . Below this separation Mi is still 
surrounded by the sea of field stars of the core. The density 
enhancement of surrounding stars in its wake after gravi- 
tational interaction further decreases the velocity of M2 so 
that it keeps on spiralling inwards due to dynamical fric- 
tion. The kinetic energy is lost to these stars, resulting in 
a heating of the core, but the stars do not gain enough en- 
ergy as to leave the system. Eventually, after some 10 7 yr, 
the distance between both the BHs shrinks to the cusp- 
radius of Mi, r cusp = r core Mi/M corc , where the velocity 
dispersion of the stars of the core equals the keplerian 



Since a BHB has a non-spherical symmetry it might velocity in the potential of Mi (Begelman, Blandford, & 


cause ; 
to poss 


,n initially spherical symmetric stellar distribution Rces 198C; 


Mikkola & Valtonen 1992 




Vecchio, Colpi, & 


ibly assume a torus-like structure. Hence, to find an Polnarev 1994 


; Quinlan 1996 


). Inside r cusp the kinematics 



answer to both the questions we performed numerical sim- 
ulations of the stars in the potential of the BHB by solving 
the three-body problem for each star of the initial distri- 
bution. Integrating over all stars of the cluster we obtain 
the properties of the stellar population and can trail its 
evolution under the influence of the binary. Following the 



suggestion by Kazanas (1989), Alexander & Netzer (1994) 



and Alexander & Netzer (1997) that the BLR is made up 



by stars and their winds, we propose that the gas and 
dust of the torus is bound to the stars within their winds. 
The dynamics of these stars in the potential of the binary 
then supports the torus in the vertical direction, solving 
the problem of how to keep the torus geometrically thick. 

In this paper we investigate the conditions for which 
the bl ack holes merge and how long such a merger lasts. 



of the stars and M2 is dominated by the potential of Mi , 
and the stars are now interacting with both BHs rather 
than with Mi only. This is the point where we start our 
calculations. Neglecting the influence of the other stars of 
the core we are now dealing with the three-body prob- 
lem for a single star in the potential of both BHs. In the 
interaction with the black hole binary (BHB) some stars 
gain enough energy in order to be ejected from the central 
region. On the other hand the more distant stars just per- 
turb the binary's center of mass and leave its semi-major 
axis unaffected and the ejection of individual stars be- 
comes the dominant process for further hardening. Thus 
stars moving on orbits with radii of about the semi-major 
axis of the binary contribute most to the shr inking of the 



binary in this stage. This is also found by Hemsendorf. 



We wi ll first use a mass-ratio of q = 10 for the black Sigurdsson, fc Spurzcm (200l| ) in their numerical calcula- 



holcs to discuss the results in more detail. Afterwards we 
compare them with those obtained for the ratios q = 1 
and 100. 

In a following paper (paper 2) we will concentrate on 
the configuration of the stars remaining bound to the black 
holes, their ability to constitute a dusty torus and the 
observational consequences in order to answer the second 
question. 

2. The model 

The evolution of the orbit of a comparatively low mass 
secondary black hole moving around a massive BH Mi in 
the center of a massive galaxy is investigated with ana- 
lytical and numerical means by Polnarev & Rees (1994). 



tions. If there are sufficient such stars the hardening con- 
tinues till eventually the black holes coalesce due to the 
emission of gravitational radiation. Otherwise the harden- 
ing stalls if the inner region is not sufficiently refilled with 
matter by some other process. 

Therefore, once the black holes approached each other 
as close as the cusp radius, the eccentricity of their orbits 



is likely to be very small (Milosavljevic & Merritt 2001) 



Dynamical friction between both galaxies results in a loss 
of energy and angular momentum so that the cores of the 
galaxies quickly approach each other. The probability of 
merging with a given initial value of pericenter is roughly 
proportional to the square of this value. Thus small initial 
pericenters are unlikely and in the most cases Polnarev fc 



so that we assume them to move on circular orbits around 
each other. 



2.1. Collisionless systems 

The net force acting on a star in a galaxy is mainly deter- 
mined by the large structure of the galaxy. Consequently 
the gravitational force of the mass distribution in the 
galaxy varies smoothly with space as it acts on a single 
star. 

The number n re iax, how often an individual star of 
mass m has to cross a system of N identical stars so that 
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the stars velocity v is changed by order of itself is 



^rcla 



N 
8 In A 



(1) 



The parameter A may be written as A = i?/& m in ~ 
Rv 2 /Gm w N, see |Binney fc Tremainc (1987| ). The largest 
possible impact parameter is limited by the scale R of the 
stellar distribution. If, on the other hand, the impact pa- 
rameter falls below the limit & m j n = Gm/v 2 the assump- 
tion of nearly a straight-line trajectory, made to obtain 
expression (|l|), breaks down. 

The time a star needs to cross the volume of the stellar 
distribution is of order t cross = R/v so that the relaxation 
time may be defined as t ro i a x = ftrelax Across- Now, for a 
total amount of N = 10 s solar mass stars we obtain for 
the number of crossing times u re iax ~ 6.8 x 10 5 . If we 
assume the star to be of solar mass and the cluster to 
extend to radii of ~ 5 pc we obtain the typical velocity of 
a star to be v = ^GNm/R w 300km/s. This allows to 
compute the relaxation time, which is t r dax ~ 10 10 yr. 

Thus, compared to the trajectory a star would follow if 
the other matter would be perfectly smoothly distributed, 
individual stellar encounters perturb this trajectory only 
over of order 0.1 AT/ IniV crossing times. This means that it 
takes several crossing times to deflect a star from its mean 
trajectory and therefore it is possible to obtain some un- 
derstanding of the dynamics by investigating the orbits of 
the stars in a suitable mean potential without taking into 
account individual stellar encounters. Since the time range 
of our simulation is smaller by a factor of ~ 10 than the 
relaxation time we can neglect individual stellar encoun- 
ters in our calculations. Thus the evolution of a stellar 
distribution in the potential of a BHB can be modelled 
by solving the equations of motion of the individual stars 
of a cluster separately. Afterwards the results of the sin- 
gle stars can be combined to model the evolution of the 
complete cluster. 



2.2. Restricted three body system 

If we assume that the potential in the central region is 
dominated by the two black holes we can neglect the mean 
potential of the stellar distribution when we are solving 
the equations of motion for the stars. Since the cluster 
can be treated as a collisionless system, we can solve the 
equations of motion for each star separately, i.e. we are 
dealing with the restricted three body problem. In the 
following all coordinates are given relative to the center of 
mass, which is identified with the origin of the coordinate 
system. The axis of rotation is the z-axis and we always 
assume Mi > M 2 and correspondingly for the mass-ratio 
q = M1/M2 > 1. The two black holes are moving on 
circular orbits around each other in a fixed distance of 
a = 1 pc. 



To write the equations of motion in a dimensionless 
form we used the following normalization constants 



ro = 
to = 

$n = 



a = 1 pc , 



V G{M\ + Af a ) 
1492 



G{M 1 + M 2 ) 



q _A 3 

M 8 2 a^ c yr, 



(2) 



J/kg, 



with the mass of the primary BH Mi being fixed to 
1O 8 M (see appendix g). M 8 = Mi/1O 8 M is the di- 
mensionless mass of the primary BH in units of 10 8 so- 
lar masses and equal to one throughout this paper. The 
quantities a pc and r pc denote the dimensionless distance 
between the BHs and the radial distance to the center of 
mass respectively, both scaled to one parsec. In the follow- 
ing we indicate dimensionless quantities with a tilde ' ~ ' 
on top. 

The numerical integration of the equations is processed 
faster in the comoving frame where the BH masses are 
stationary on the cc-axis and the rotation axis is pointing 
along the z-axis with the z = 0-plane being the equatorial 
plane of the BHs. In this system the equations of motion 
read (see appendix |X]): 



x 



V x = - 



X — X\ 

l + q \ q \r - f 1 1 3 



1 + q \ q \r - fi| 3 



1 + q \ \r - n 



X + qx\ 
|r + qfi| 3 



I? + <?ri| 3 



I? + <?fi| 3 



+ 2v y +x (3) 
-2v x +y 



These have been solved numerically after being supplied 
with a set of initial conditions for the stars which we will 
discuss in the next section. 

In order to check the accuracy of the numerical inte- 
gration, we keep track of the deviations of the Jacobian 
Integral (appendix |b|), which is a conserved quantity in 
this problem. 

2.3. Initial stellar distribution in phase-space 

The surface-brightness profiles of many elliptical galaxies 
are well fitted by an isothermal sphere out to a few core 
radii. On the other hand rotation curves of spiral galaxies 
are often remarkably flat out to great radii, and this sug- 
gests that we should construct models that deviate from 
the isothermal sphere only far from their cores (Mihalas 
|fc Binney 198l| ; [Binncy fc Tremainc 1987| ). 



The structure of an isothermal self-gravitating sphere 
of gas is identical with that of a collisionless system of stars 
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whose density in phase-space is given by Eq. (C.4), see ap- 
pendix |^. This correspondence between the gaseous and 
stellar-dynamical isothermal sphere originates in the ve- 
locity distribution of both systems. Integrating Eq. ( |C.4[ ) 
over the volume yields the Maxwell-Boltzmann distribu- 
tion for the velocities. As we know from kinetic gas the- 
ory this is the equilibrium distribution for a gas whose 
particles are allowed to collide elastically with each other. 
Therefore, if a system of particles is described by a distri- 



Number density of stars p(f) 



bution function according to Eq. (C.4), it makes no dif- 
ference whether the particles collide with each other or 
not. This makes the isothermal sphere a very well suited 
initial distribution for our purposes. The mean and the 
mean-square velocity of the stars in the isothermal sphere, 
(v) = y^S/na and (v 2 ) = 3er 2 , are independent of the posi- 
tion. Since the velocity dispersion is isotropic everywhere 



its components are the same: 



Interestingly Milosavljevic & Merritt (2001) 



= <T . 

find in 

their numerical calculations that the radial density pro- 
files of the pre-merger galaxies and of the merged galaxy 
just after the formation of a hard BHB are essentially the 
same for a short time. For the initial density distribution 
they used a powerlaw with index —2. 

In our model we distribute the star s according to 
the singular isothermal sphere (Eq. (C.9)) in the range 
0.1 < f < 50. In order to obtain the Maxwell-Boltzmann 
distribution in velocity each component has been dis- 
tributed according to a gaussian. Since we are interested 
in stars which initially are bound by the potential of the 
BHB (E < 0) and do not leave the system immediately 
after the calculations started we choose the free parame- 
ter a, the velocity dispersion, to be a third of the escape 
velocity, a = v csc /3 = V2$/3. 

The initial conditions for 8000 stars have been gener- 
ated and supplied to the Eqs. (|^). These are solved with 
a Runge-Kutta code of fourth order based on the routine 
rkck by Press et al. (1995) ) . The stepsize is selfadapting to 
ensure that the desired accuracy is always maintained and 
that the stepsize does not become too small in order to 
save computing time. Due to our choice of the initial con- 
ditions all stars are bound to the binary in the beginning, 
having negative energies. Some of them are ejected later 
during the run time. We define a star as being ejected if 
the following three conditions are fulfilled simultaneously: 

— the energy of the star is positive (E > 0) 

— its radial distance in orbit is bigger than 50 a 

— the radial velocity component is positive (v r > 0). 

The simulation is stopped after the time i max = 
5 x 10 5 (in dimcnsionlcss units) is reached or if the 
star has been ejected before. i max corresponds to 
about 80000 re volutions of the BHB or to w 7.46 x 
10 8 M- 1/2 4/ c V<7/(<Z + l)yr. 

In the next section we present the results we ob- 
tained for the simulation using a mass-ratio q = 10 
(M 2 = 1O 7 M ). Afterwards we will compare them with 
the results we obtained for q = 1 and q = 100. 



0.01 + 



0.0001 t 



le-06 t 



le-08 i 



le-10 f 



TP, t = 
EP, t = 
BP, 1=0 

BP, t — ^max 



le-12 




0.1 



10 



100 



1000 



Fig. 1. Initially the central region is dominated by the 
ejected population (EP), the outer regions by the bound 
population (BP). Finally, at i max , the bound star distri- 
bution follows a powerlaw with index —4 in the heated 
region (f > 50), while an index ~ —2 is maintained in 
the range 10 < f < 50. The inner parts are scoured out 
and a maximum emerges at f ~ 3, showing the torus-like 
configuration the stars assume. For f < 2 a cusp of stars 
bound to Mi only is left, f is normalized to 1 pc and p 
such that the area underneath the solid line is 1 . 

3. Results for the singular isothermal sphere 

In the following all quantities are given in and related to 
the observers frame if not indicated otherwise. The origin 
of the coordinate system is the center of mass of the BHB, 
with its rotation axis pointing along the z-axis. At t = 
both black holes lie on the x-axis, with Mi situated on the 
positive and Mi on the negative seminaxis. The y-axis is 
perpendicular to both the others so that all axes form a 
right handed tripod. In spherical coordinates 9 denotes 
the angle to the positive z-axis, and (j) is the azimuth in 
the x — y plane. 

A key feature is the torque exerted by the rotating 
black holes on the stars, which changes the star's orbit 
and can eject it from the system. 

During the simulation we kept track of the single stars 
what enables us at all times to assign them either to the 
ejected population (EP), which is constituted by stars be- 
ing ejected before the end of the simulation, or to the 
bound population (BP) when they remain tied to the bi- 
nary. Both these populations combine to give the 'total 
population' (TP). Quantities referring to EP, BP and TP 
are denoted with the indices 'ej', 'bn' and 'total' respec- 
tively. 




x=0 plane z=0 plane 

Fig. 2. Cuts through the stellar density in the comoving frame are shown with contours logarithmically scaled. The 
right panel displays the equatorial plane (BHs marked by black spots). Perpendicular to it the x — 0-plane is shown 
(left panel), with the y-axis drawn as dashed line, so that the BHs are in front and behind the paper-plane. The initial 
distribution is a gaussian and the massratio is 10. Time increases from top to bottom (indicated by t). After initially 
stars close to M2's orbit (w 1 pc) are ejected (2nd row) also the polar regions are depleted and a torus in the equatorial 
plane finally emerges at r ~ 3pc (3rd row). 
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3.1. The dynamics of the stars 

The density distribution of both populations EP and BP 
shows that the bound stars dominate the outer regions 
^ 10, dashed-dotted line in Fig. |l|). With decreasing 
radius the fraction of ejected stars (long-dashed line) is in- 
creasing with both populations being comparable in num- 
ber density in a distance f ~ 10, corresponding to 10 pc. 
At f w 3 the density of the BP has a maximum and after 
a gap in the distribution between 1 and 2 pc, where M2 is 
circling around Mi , the density of bound stars is increas- 



Torque divided by radius, f 
I 



1^ 



ing again towards smaller radii. Niemeyer & Biermann 
1993| ) could reproduce the range of observed mm to in- 



frared spectra of radio-weak quasars. In their model dust 
is confined to molecular clouds in a disk configuration and 
is heated by the central engine. According to the model by 
Barvainis (1987 ) the dust is heated to its evaporation tem- 
perature by the central UV luminosity in ~ 1 pc distance. 
This eva poration radius is taken as the inner radius of the 
torus by Lawrence (199l|) and is in very good agreement 
with the inner edge we find for the bound stellar distri- 
bution at about 2pc (Fig. |). |Krolik fc Bcgclman (1988|) 
determine the inner radius of the torus from the balance 
between cloud evaporation by central radiation and in- 
flow by dissipative processes in the torus and also obtain 
~ 1 pc. At distances smaller than one parsec (f < 1) only 
62 stars of the BP (a fraction of 1/80) are found at t = 0. 
These are bound to the primary black hole only, and till 
the end of the simulation the number of stars in this cusp- 
region is increased by just one which has been captured 
by Ml 

The evolution of the initial density distribution (a 
gaussian) to its final state is shown in Fig. |2[ It displays 
cuts through the 3-dimensional density distribution in the 
comoving frame of the binary. The right panel shows the 
equatorial plane (z = 0) with the BHs marked by the black 
points with white frames. Ml to the right. The cuts in the 
left pahel arc perpendicular to the equatorial plane and 
contain the x = 0-plane with the y-axis indicated by the 
dashed line. The time proceeds from top to bottom and 
corresponds to 10 4 86 , 10 657 and 10 7 71 yr respectively for 
(7=10. The first row displays the distribution very close to 
the initial state. After 10 6 57 yr (second row) we see that 
mainly stars close to the secondary's orbit (f — 1) have 
been ejected and a gap appears at this distance in the 
equatorial plane leaving a shell in the range 1.5 < f < 3 
behind. At about f — 2 the density has a maximum and 
we can detect a shallow torus emerging at this radius. 
But still the polar regions arc populated as dense as the 
equatorial plane, see Fig. |[ second row. 

However, as is illustrated in the 3rd row of this figure, 
after about 10 7 71 yr the central regions have been scoured 
out and the hole in the distribution has a cylindrical ge- 
ometry elongated along the binaries rotation axis. In the 
equatorial plane finally a torus with a radius of about 3 pc 
emerges, showing that the density distribution develops 
from the initial sphere over a shell-like distribution (2nd 
row) to the final torus. Thus one important condition of 




Fig. 3. The component L Zt x/f = f(j)sm 2 (9) of the nor- 
malized torque of the binary, which acts on a star moving 
on a fixed circular orbit with radius f is displayed as func- 
tion of time. The angle enclosed by the rotation-axis of the 
binary and the orbital spin-axis of the star is denoted by 
9 and is varied in the range from 0° to 90°. Perturbations 
are stronger for orbits through the polar regions; these 
stronger perturbations lead to secular changes in the or- 
bit, finally leading to ejection for about half of the stars 
in the polar cap (see Fig. ||). 



the unification scheme is fulfilled, namely that the stars 
do assume a torus-like distribution at a distance which is 
in agree ment with the or i gin of thermal infrared emission 
by dust (Barvainis 1987 ; Krolik fc Bcgclman 1988 ; Haas 



et al. 2000). Whether this torus also satisfies the other re- 



quirements, like optical thickness, will be investigated in 
our paper 2. 

The basic topology of the final density distribution 
can be understood in physical terms, where we consider 
the torque exerted by the two black holes on an orbiting 
star. Fig. || shows of the normalized torque the component 

L z ,i/r — r<j) sin 2 (9) relative to the z-axis of the BHB as 
a function of time for different angles 9 between the rota- 
tion axis of the BHB and the symmetry axis of the star's 
orbit. The bigger the angle 9 is (i.e. for orbits through 
the polar regions), the more is the star's trajectory dis- 
turbed by the influence of the two BHs. The curves are 
symmetric since for simplicity it has been assumed that 
the star is moving on a keplerian circle of radius R around 
a point mass which corresponds to that of the binary and 
is located at the origin. Of course real trajectories are dis- 
turbed by the torque and the curves will not be symmetric 
any more. The cumulative effect of these large excursions 
in the polarcap regions deplete the stellar population leav- 
ing a torus behind. 
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Fig. 4. The angular momentum clearly separates the BP 
(black dots) from the EP (grey dots), which initially have 
less angular momentum. The line of transition from the 
EP to the BP is drawn by eye. Its almost constant value 
is decreasing with increasing mass-ratio q. Along this line 
both populations diffuse into each other, showing that the 
transition is smooth. 



Because of the conservation of the Jacobian Integral 
(Eq. (B.4)) the stars which gain energy also gain angular 
momentum in the z-componcnt and thus help to enhance 
the density in the equatorial plane. In paper 2 we will 
discuss the distribution of the BP and investigate in the 
properties of the torus. 

At t = the finally ejected stars can be clearly distin- 
guished from the bound population by their total angular 
momentum L or pericenter f _ . The mean value of the an- 
gular momentum of the EP is smaller by a factor of 3.22 
compared to the BP (see Table lion page 19). This separa- 
tion is clearly illustrated in Fig. |4l where we have plotted 
the initial values of angular momentum of the stars versus 
their radius. The ejected stars, marked by grey dots, have 
low angular momenta, while stars staying bound to the 
binary till the end of the simulation (black dots) exceed 
a certain value. Both populations can be separated from 
each other by a transition line, as seen in Fig. ^. The angu- 
lar momentum has been normalized to the value of a star 
which moves on circular orbits in a distance corresponding 
to the separation of the BHs (a = 1 pc) around the com- 
bined mass of the BHs Mi + Mi which is situated at the 
center (L — ma 2 /t = m^jGMxa{l + q)/q)- Along this 



transition-line, which is drawn by eye and approaches a 
value of about 2 in these dimensionless units, the different 
populations diffuse into each other so that the transition 
from one population to the other is smooth. 

This value is double the angular momentum of a star 
moving on a circular orbit with the radius corresponding 
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Fig. 5. Ejected stars (grey dots) cover a region with f_ < 
2 and f+ > 0.7, marked by the horizontal and vertical 
lines. This allows them to closely approach the orbit of 
M 2 . Bound stars (black dots) avoid this unstable region. 
One unit corresponds to 1 pc. 



to the separation of the binary. Or, the other way round, 
this transition-value of L corresponds to an orbit with 
a radius four times as big as the semi-major axis of the 
binary. 

However, due to their much lower angular momenta 
the finally ejected stars come much closer to the black 
holes. For the mean pericenters we get (r C j_) = 0.86 and 
(rbn-) = 10.42 showing that they differ by more than a 
factor of 10. In Fig. [5] we have plotted the initial apoc- 
enters on the cc-axis versus the initial pericenters on the 
y-axis. Without loss of information and to keep the figure 
distinct, we plotted only 1/4 of each population, chosen 
by random. While the pericenters of the ejected stars stay 
below ~ 2 (i.e. 2pc, horizontal line), their apocenters ex- 
ceed ~ 0.7 (0.7 pc, vertical line). Hence they cover a region 
which always allows them to come very close to the orbit 
of the secondary BH, which rotates at a distance of 0.9 pc 
around the center of mass. Here the stars undergo violent 
interactions with the binary and eventually gain sufficient 
energy (E(t max ) > 0) to be kicked out. On the other hand 
the bound stars avoid this region and stay at distances 
large enough to avoid such violent interactions with the 
binary. Only very few stars, less than 1.3% from the BP, 
have apocenters less than ~ 0.7 pc. These are bound to the 
primary black hole Mi only and are not ejected because 
the influence of the secondary BH is too small. They form 
the small density peak around Mi which is seen in Fig. |[ 

For the same apocenter the ejected stars have on av- 
erage smaller pericenters than the bound stars and so are 
moving on orbits with higher eccentricities. As a conse- 
quence the radial velocity component should exceed the 
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Fig. 6. The BP is shown to be strongly tangential 
anisotropic while the EP is radially anisotropic. As the 
eccentricity the velocity anisotropy (3 increases with f for 
all populations but for the TP which is initially set to zero 
at all distances. At t max the BP is radially anisotropic in 
the heated region f > 50. In the cusp region f < 2 no clear 
tendency is detectable. 

tangential component, leading to a radially anisotropic 
velocity. But due to our choice of initial conditions the 
velocity anisotropy 



[3 = 1 



(vl) 



(4) 



is initially zero at all radii for the total population, as can 
be seen in Fig. [| (solid line) . The EP is radially anisotropic 
for all radii (/3 > 0, long-dashed line) and the velocity 
anisotropy is increasing with distance. While (3 is also 
increasing with radius for the BP, it is always tangen- 
tially anisotropic ((3 < 0, dashed-dotted line), the more 



the smaller the radius is. Quinlan & Hcrnquist (1997) also 
find (3 of the bound stars to grow with f and obtain in 
their numerical experiments a minimum value of — 1 , also 
starting with (3 = at all radii. They suspect the velocity 
anisotropy to decrease further towards smaller radii. This 
is what we find in Fig. ^ where [3 becomes as small as 
—2.5 at a distance f ~4 (—3 for q — 1), the dotted line. 
This is much smaller than the values —0.3, predicted by 
the adiabatic-growth model for the formation of massive 
black holes ( Quinlan, Hcrnquist, fc Sigurdsson 1995|), and 



—0.4 obtained in numerical simulations by Milosavljevic 
& Mer j-itt (2001]) . 

eccentricity of 



The 



both 

-)/(?+ 



populations e = 
f- r_) is increasing 



VI + 2EL 2 = (r+ - r 
with distance from the binary. Close to the center the 
bound stars have to move on almost circular orbits in 
order not to come closer at the pericenter and to avoid 



strong interactions with the black holes. With increasing 
distance also the eccentricity can increase as long as 
a big enough pericenter (f _ > 2 pc) is maintained i.e. 
not smaller than twice the separation of the BHs. The 
EP must have a pericenter below this value so that its 
stars can interact violently enough with the binary at 
the point of the closest approch to the black holes in 
order to become ejected. Thus with increasing apocenters 
also their eccentricity increases, where e — * 1 as f — > oo. 
Consequently these stars have low angular momenta 
and their orbits become the more radially anisotropic 
the bigger their distance to the center is. Hence their 
kinetic energy is stored in radial rather than in tangential 
motion. 



3.2. Loss of L and subsequent merging of the black 
holes 

As time proceeds there are almost no changes in the mean 
value of the total angular momentum (L) for both popu- 
lations, BP and EP. With the angular momentum being 
normalized to Lq — ma 2 /to, where m = M Q is the mass of 
a star, the mean values of L z are —0.045 and 0.063 for the 
BP and EP respectively. While on average the bound stars 
are slightly counterrotating, the ejected stars are corotat- 
ing so that the total population shows no net rotation 
in the beginning. This matches the results obtained by 
[nnancn (1979| ) and |lnnanen (1980| ). The mean angular 
momentum (L z ) does not change for the bound stars as 
time proceeds, but it increases by a factor of about 4.54 
for the ejected population (see Table pi page 19). Such 



a growth is expected, because of the conservation of the 
Jacobian Integral (Eq. (B.4)). All ejected stars gain en- 
ergy and consequently their angular momentum in the 
z-component has to be increased, at the expense of the 
BHB. Hence the ejection of stars leads to a net loss of 
angular momentum and thus to hardening of the binary. 

In order to get some quantitative information about 
the hardening we can use the statistical values tabulated 
in Table [l]. They allow to compute the mean angular mo- 
mentum that is extracted from the binary by a single star 
and thus also the number of stars which is required to 
carry away all the binaries angular momentum so that 
the black holes coalesce. But, since we made certain as- 
sumptions and neglected some effects, these numbers will 
give just an order of magnitude estimate. On the one hand 
the feedback of the stars on the binary has not been taken 
into account because we fixed the distance of the black 
holes throughout all calculations. As a consequence the 
"cross-section" for the stars to be ejected is not shrinking 
with increasing time, when ejected stars are extracting 
angular momentum from the binary. Thus the fraction of 
ejected stars, obtained from the simulation, is bigger and 
the initial total number of stars needed to allow for merg- 
ing will be a lower limit. For the same reasons also the 
time-scale of the merger of the BHs, which we will esti- 
mate in this section, will be shorter than in reality. On 
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the other hand we neglected star-star interactions, which 
would support the loss-cone feeding and therefore increase 
the fraction of ejected stars. Also we assumed the black 
holes to move on circular orbits. If the orbits were ec- 
centric, gravitational radiation would become important 
earlier at a bigger semi-major axis resulting in a faster 
merger and less stars needed in the cluster. Both the lat- 
ter effects counteract neglection of the shrinking of the 
black holes' distance, but as shown previously they are of 
minor importance. ' 
However, the angular momentum of the BHB is 



Angular momentum loss rate =kf-, q = 10 



L„ = 10 8 M f 

= 10 8 M 8 (l + g)(fi x h) 



1 



ri x n + -r 2 x r 2 

q 



10 8 M s a 



1/2. 



1 



1 



<1 



AL Z 



1 10 8 Af 8 
1 + q AL Z 



4.07 x 10 7 M 8 



N h 



0.001 cr 



0.0001 - 



le-05 - 



le-06 - 



(5) 



where we used the relation f 2 = —qr± between the 
positions of the two black holes. With the expressions 
|?i x f 1 1 = |fi||fi| and 5 = |f x - r 2 | = (1 + g)|?i| = 1 the 
absolute value of the angular momentum becomes 



(6) 



If AL Z is the average angular momentum extracted from 
the BHB by a single ejected star, the number of ejected 
stars needed to carry away all the binary's angular mo- 
mentum amounts to 



(7) 



In the last step we used the values q — 10 and AL Z = 
(L z (t max )) - (L z (t = 0)) « 0.22 taken from Table @. With 
the ratio of bound and ejected stars N-t, n /N e j = 1.63 from 
Table [l] the required total number of stars to let the BHs 
merge is: 



1.07 x 10 8 . (8) 



If, as we assume, every star has on average the mass of 
the sun, the initial star-cluster is as massive as the binary 
(Mi + M 2 = Mi(l + 1/q) = 1.1 x 10 s M Q ). Hence our 
neglection of the mean potential, generated by the stellar 
cluster compared to the potential of the BHB close to the 
center, is justified. 

To get an estimate of the timescales on which the black 
holes merge and which serve as a lower limit as stated pre- 
viously, we have to know about the loss-rate of the angu- 
lar momentum of the binary. This is displayed in Fig. [?] as 
function of time. The data are normalized to the number 
of simulated stars so that the curve represents the mean 
rate at which a total of one solar mass extracts angular 
momentum from the binary when ejected gradually till 
the end of the simulation. This mass will be referred to in 
the following as a "representative star" . A simple fit to the 
data enables us to estimate the time needed for the black 
holes to merge and also to calculate the number of stars 
required so that the black holes can coalesce. This is then 




le+06 



Fig. 7. The loss-rate of the binary's angular momentum 
is largest in the beginning. For t > 1000 the loss-rate 
approximates a powerlaw with index w —3/2, shown by 
the dotted line (fit 2). The dashed curve (fit 1) is used 
to calculate the number of ejected stars. One time unit 
corresponds to 1423 yr. 



compared with the number iV C j obtained above by statis- 
tical means and should be the same if the fit is sufficiently 
accurate. 

After a delay of t w 26 the loss-rate is increasing 
steeply, passing through a maximum at t ~ 150 (corre- 
sponding to 2.1 x 10 5 yr) and then quickly approaches a 
powerlaw with index ~ —3/2. This initial delay is caused 
by two reasons: first the stars have to interact with the 
binary before they gain enough energy to be kicked out 
since all stars are bound at the beginning. Afterwards the 
star, having now a positive energy, has to move to a dis- 
tance of 50 pc in order to be registered as ejected by the 
code. The exponential decrease of the loss-rate for t > 10 3 
(~ 1.4 x 10 6 yr) shows that the system evolves faster in the 
early stages and the angular momentum varies very slowly 
at later times. If the black holes get sufficiently close due 
to the ejection of stars, gravitational radiation will be- 
come strong enough to complete the merger. For t > 10 4 
the curve is not very smooth because of the small number 
statistics. As we will see later this corresponds to the time 
required by the BHs in order to merge, ~ 1.4 x 10 7 yr. The 
lossrate can be approximated by a powerlaw with index 
—a — —3/2, multiplied with a function which cuts it off 
at t d = 20: 



dL 



dt 



£ = ki- 



rn 
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Integration with respect to time yields 



C(i) = 



dL, 



-+dt = 



dt 




where F is the hypergeometric function. As time tends 
to infinity C(t) tends to zero. The angular momentum 
extracted from the BHB by such a representative star in 
its dependence on time is C(t) — £,(t&). The integration of 
Eq. (^) starts at ?<j to take into account the delay till a star 
is registered as ejected. The average amount of angular 
momentum the two black holes can lose to the mass of 
the representative star is 



C(t -* oo) - C(t d ) = -£{t d ) 

2k r(i + ^)r(i + «) 



r(i 



(ii) 



^ + «) 

where T denotes the Gamma function. To compute the 
fit-parameters we use an implementation of the non-linear 
least-squares Marquardt-Levenberg algorithm and obtain 



1.77, 77 = 1.18 and 



9.84. 



(12) 



These parameters provide a sufficiently good approxima- 
tion to the data as can be seen in Fig. |?], where the fit 
is displayed by the dashed curve (fit 1). Substituting the 
parameter values to Eq. ([ll]) results in 

(13) 



Loo = 0.258 . 



This is about a fourth of the angular momentum of a star 
circling at 1 pc distance around a point mass correspond- 
ing to that of the binary. Now the number of ejected stars 
required to carry away all the binary's angular momentum 
is obtained in an analoguous way as in Eq. (Q): 



N f .. 



L. 



1 10 8 M 8 
1 I 9 Loo 



(14) 



Comparison of both expressions shows that the average 
angular momentum gained by a single star, AL Z , has been 
replaced by the angular momentum which is extracted on 
average by one solar mass till the end of the simulation, 
Loo. For the same ratio Nb n /N e j as above we now obtain 
the numbers 



iVoj « 3.47 x 10' 



and JVtotai w 9.33 x 10 7 



(15) 



They are in very good agreement with the numbers we 
got above by statistical arguments only, confirming our 
finding that a comparable amount of mass in stars as in the 
black holes is required to remove all the binary's angular 
momentum. 

Comparing both the loss-rates of angular momentum 
due to gravitational radiation and ejection of stars allows 
us to compute the time and distance where gravitational 
radiation eventually becomes the dominant physical pro- 
cess for the merging of the black holes. As before we as- 
sume the number of stars to be sufficient so that all an- 
gular momentum of the binary is lost as time tends to 



infinity. Hence the angular momentum lost by the binary 
as a function of time is 



ii«>t(£) = N ej (C(t) - C(t d )) = ~^-(£(i) + Loo) 



Ufa) + i 



\L C 



(16) 



For t < t d we assume all quantities to be constant. The 
angular momentum left to the two black holes is simply 
its initial value diminished by Ljost j an d because of its de- 
pendency on the distance a(t) (Eq. (0)) we get the relation 



L m (t) = L m {i A )y 



L m (i&) — Liost(i) 



-Cii) 



(17) 



The separation of the black holes which we kept fixed to 
1 pc (a = 1) during the simulation, is now treated as func- 
tion of time with an initial value d(t < t d ) — 1. Solving 
Eq. ( |l7| ) for the separation, a(t) — (C(t) /Loo) 2 , the shrink- 
ing rate of the binary due to ejection of stars can be writ- 
ten as 



da _ £{t) dC(t) 
dt Loo dt 



(18) 



This rate has now to be compared with the shrink- 
ing rate caused by the emission of gravitational radiation. 
Using Mtotai = Mi(l + q)/q and fj, = M 1 /(l + q), the total 
and reduced mass respectively, the energy of the binary 
reads 



E m 



1 GAtM tota i l 
E 2r a 



1 



(19) 



E m being dimensionless and E = GfiM tota \/r the nor- 
malization constant. The en ergy loss of the binary due 
to gr avitational radiation is ( Misncr, Thome, fc Wheeler 



1973|) 



dE m 

dt 



32 G[i 2 rfo 5 
~5 E Q c 5 
32 / GMi 



q+l 



(20) 



Forming the time derivative of Eq. ( {l9| ) and with the help 
of Eq. (|20|) the shrinking rate of the binary results in 



da 
dt 



64 
5c 5 

1 



GMi 

q + l 

Q 3 



q + i 



--3 



S.3 



(21) 



For convenience we introduced the dimensionless 'shrink- 
ing time' t s h in the last step. Multiplied with to this ref- 
erence time scales as 



t Bh t = 2.33 x 10 ib M, 



-3 4 
pc 



yr 



(22) 
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Fig. 8. The shrinking rates due to ejection of stars (solid 
line) and emission of gravitational radiation (dashed line) 
are shown to have very different slopes. Therefore the re- 
gion of the transition is quite narrow. The separation a is 
normalized to 1 pc. 

Integrating Eq. ( |2l| ) from t = 0, when the distance is ao, 
to i yields 



a(t) 

a 



1 - 



sh 



(23) 



We define a tr as the distance of the transition where emis- 
sion of gravitational radiation becomes more effective than 
ejection of stars in extracting angular momentum from 
the binary. If we now assume that the integration of the 
shrinking rate ( ^l| ) starts when the BHs happen to be just 
this distance apart (a = a tr ), we get 



5(f) = A _J_ 

a tr V £mt 



In this expression 



t — _il 



1 



(24) 



(25) 



is the time the black holes still need to merge completely 
once gravitational radiation starts at a = 2it r to govern 
the merging process. 

If an in Eq. ( p3| ) is equal to 1, i.e. the separation where 
we started the simulation and where the BHs are 1 pc 
apart from each other, the time needed to reduce a by 
some large fraction because of emission of gravitational 
radiation only is of order of t sri in ~ 1-56 x 10 12 to w 
2 x 10 15 yr. Thus in the beginning, at distances as large 
as 1 pc, the hardening of the binary caused by emission of 



gravitational radiation is completely negligible and pro- 
ceeds on time-scales orders of magnitudes higher than 
shrinking due to ejection of stars. To compare these two 
rates of hardening we have to rewrite Eq. (^) in its depen- 
dence on distance a instead of t. Unfortunately, because 
of its complicated form, we can not build the inverse func- 
tion of a{t). However, to proceed we first have to further 
simplify our fit function. This is done by resetting the pa- 
rameter k to zero, so that the cut-off funtion in Eq. (JsJ) 
vanishes and we are left with a pure power law only. In 
order not to have too large differences between this fit and 
the data for small times and to take account of the initial 
increase we shift the parameter from 20 to 50. Hence 
the system function ( |io| ) becomes 



C(t) 



71 — a 



(26) 



To maintain the number of ejected stars needed for a final 
merger of the black holes we recalculate the parameter k 
using Eq. (|ll|). Thus = kth~ a /(a — 1) and therefore 



we get k — L 



•i)*r 



0.91. Since we are just inter- 



ested in an estimate of the timescales of the merger, these 
parameters still provide a sufficiently good fit to the data 
as can be seen in Fig. [?], where the power law is shown 
as dotted line an denoted as 'fit 2'. Now the distance as 
function of time, a = (td/t) 2 ^ a ^ 1 \ can be solved for i. 
With a — 3/2 the shrinking-ratc da^ 1 /dt is constant, as 
is found by Milosavljevic fc Merritt (200l| ) in N-body cal- 
culations. Substituting the time as function of the separa- 
tion together with Eq. (|2^) in the the shrinking rate due 
to star ejection (fLST) yields 



da 
dt 



2(a-l) = 



a 2( Q -l) 



(27) 



Equating this expression with the hardening rate caused 
by emission of gravitational waves ( pl| ) we can solve for the 
distance of the transition ct tr where gravitational radiation 
takes over the process of shrinking, 



Otr 



2(a- 1) V Q 



3 + 1 



' id 



t. 



.sh 



(28) 



The merging for a > a tr is dominated by the ejection 
of stars. Since both shrinking rates have very different 
slopes (see Fig. ||) dn does not depend sensitively on the 
parameters we chose for the fit and the transition occurs 
quite suddenly over a small range in a. Substituting atr 
to a = (?d/?) 2( - a_1 ^ we can solve it for the time ftr needed 
till the distance has shrunk from its initial value a = 1 
to a = atr due to ejection of stars. For the parameters we 
used, a = 3/2 and id — 50, we finally obtain for q = 10 



1 _ , 
atr ~ Y97 = 5.1 x 10 J pc, 

t tr » 9855 = 1.40 x 10 7 yr. 



(29) 



Hence after about 10 7 yr the distance of the two black 
holes has shrunk to ~ 1/200 of its initial value lpc. 
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Merritt (200C) obtained for a tr about 1/100 of the distance 
where the binary becomes hard, which translates in phys- 
ical units to about 1/100 pc according to Milosavljcvic fc| 
Merrit ; (2001). This is in good agreement with the dis- 
tance of transition for the mass-ratio q = 1, sec Table |. 

Only the fraction ~ 1/14 of the initial angular momen- 
tum of the binary is left at this distance. With a = 3/2 
for the powerlaw index the relation between the distance 
and time of the transition is ctt r = td/ttr, and Eq. 
simplifies to 



atr 



9 + 1 

<Z 3 



/si I 



Thus the ratio of the time ranges when star ejection and 
emission of gravitational waves dominate the hardening is 
independent of the massratio: 



tt 



t 



(30) 



nit; 



For the time when gravitational radiation dominates we 
then have 



mg 



2464 = 3.51 x 10 b 



(31) 



showing that the merging process is dominated most of the 
time by ejection of stars. Quinlan (1996| ) obtains for the 
timescales needed for merging once the binary has become 
hard (a w r C usp), about 4 x 10 7 yr. This is very similar to 
our result, which is supposed to be a lower limit since in 
the simulation the shrinking of the distance between the 
black holes has not been taken into account. This time- 
scale on which the BHs merge is of the same order as the 
one we obtained for the torus to emerge, see Fig. ^. Thus, 
if a binary with a mass-ratio 10 merges due to ejection of 
stars from a surrounding cluster, a torus-shaped distribu- 
tion forms on a similar time-scale. 

The energy of the binary scales with its separation as 
5 _1 . If the distances of the stars in the cluster are in- 
creased or decreased by the same factor x as the binary's 
separation, the energy of a star also scales with this factor 
to the power of —1, i.e. E m and £J* oc x _1 . Thus the basic 
results will be unchanged as long as the binary remains 
"hard" , but the time-scales of the merging due to ejection 
of stars scale as x 3/<2 - 

The final merger of two supermassive BHs should lead 
to one of the most energetic events in the universe. As a 
result of such a merger, the spin and its direction are ex- 
pected to change (Wilson & Colbert 199!:). Consequently 
also the direction of the jet, which is aligned with the BH's 
spin, will jump into the new direction of the spin of the 
merged BH. This is a possible explanation for the observed 
X-shaped radio galaxies. This class of objects exhibits 
four jets, one pair of which shows "young" synchrotron 
emission and th e other one "old" emission, such as seen 
in B2 0828+23 flParma, Ekers, fc Fanti 1985|; |Rottmann 



Dennei -Thorpe, fc Klein 1998| ). The spectral aging of the 
secondary lobes in this object suggests that a merger has 




Fig. 9. This sketch illustrates the meaning of the different 
angles. The equatorial plane is denned by the black holes 
rotating around each other, with the angular frequency 
vector G)„ pointing along the z-axis. Mi is displayed a lit- 
tle off-center since the coordinate system is centered in the 
center of mass frame. Tilted to this plane with an angle 
#piane is the orbital plane of a star, whose angular mo- 
mentum L* encloses the angles 9l with u) m (the positive 
z-axis) and 9l_ with the negative z-axis. The line from 
the center to the current position of the star has the angle 

#nns With bJmm. 



happened about 6 x 10 7 yr ago (priv. communication with 
H. Rottmann). This timescale is very close to the one we 
have found here for the merger of two supermassive BHs 
and thus supports the idea of mergers inducing a jump of 
the jet. 

3.3. Final distribution and energy of ejected stars 



As we have shown in Sect. 3.1, both populations, the 
ejected and bound stars, can be clearly distinguished by 
their angular momentum or their pericenter. To learn 
more about their angle distribution is the task of this 
section. For this purpose we introduce some angles, as is 
illustrated in Fig. ||. We define 9l as the angle between 
the positive z-axis (the rotation axis Q m ) and the angular 
momentum of a star. The angle between L* and the 
negative z-axis is denoted by 9l_ , so that 9l_ = 180° — 8l- 
To the angle between the plane in which the star rotates 
and the equatorial plane we refer as 6* p i anc . 

The angle which is enclosed by the line connecting the 
origin with the current position of the star and the z-axis 
is denoted by 9 pos . The mean and the dispersion of this an- 
gle are very similar at t = for all three populations (EP, 
BP, TP), namely (9 pos ) » tt/2 and A9 pos « ^ 2 - 8)/4. 
These are just the values we expect for a spherically sym- 
metric distribution, as the total population at t = 0. Thus 
both populations, EP and BP, are also spherically sym- 
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Fig. 10. The density of the ejected stars at t = (spher- 
ically symmetric distributed) and after their ejection (t C j) 
as function of their polar angle is displayed. After being 
ejected the cluster is more concentrated to the equatorial 
plane 6> pos = 90°. 



metric distributed in the beginning and none of them 
shows any preferences as to certain polar angles. However, 
from Table [l] we see that the mean angular momentum (L) 
of the ejected stars increased by just a factor of 1.1, while 
(L z ) is by more than a factor 4.5 bigger at the end of 
the simulation. Thus we expect the geometry of the dis- 
tribution of the ejected population to be flattened and not 
spherically symmetric at the end of the simulation. 

In Fig. Il0| we plotted the density of the stars in de- 



pendence of the cosine of 



As stated above we see 



concentrated towards the equatorial plane of the binary. 
The same conclusion is drawn from Fig. [ll]. It displays the 
density as a function of the orientation of the orbital an- 
gular momentum of the stars, cos(#l) = L z /L. At t = 
(dashed line) the density is almost constant in cos(#l), 
thus there is no special orientation prefered for the or- 
bital plane of the stars. The density is slightly increasing 
with cos(#l), showing that the ejected stars initially are 
weakly corotating on average. After the ejection this curve 
is steeply increasing, with only few stars having L z < 
(counter-rotating) and the distribution having the maxi- 
mum at cos(#l) = L z /L — 1. This means that the an- 
gular momentum from most stars points along the rota- 
tion axis and thus their orbital plane is aligned with the 
equatorial plane. Therefore the density of the stars has a 
maximum in the equatorial plane and a minimum at the 
poles, as seen before in Fig. [Io[ The angle 6* p i ano between 



0.01 



i — i — i 1 1 r 
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within the fluctuation of the data a uniform distribution 
for t = (dashed line), according to a spherical symme- 
try. But after the stars have been ejected the density has 
a maximum at cos(# pos ) = (solid line). Thus they are ^pos 



Fig. 11. The orientation of the angular momenta of the 
EP at t = and after the ejection at i e j is displayed. They 
show a strong concentration towards large positive values, 
so that the orbital plane is close to the equatorial plane. 

orbital and equatorial plane (see Fig. ^J) is the same for 
clockwise (L z < 0) and counter-clockwise {L z > 0) rota- 
tion and is equal to the minimum of 9l and 6l_. Gaining 
now angular momentum in the z-component leads to a de- 
creasing angle 9l and increasing . A formerly counter- 
rotating star with cos(^l) = — 1 could now cross the poles 
(cos(6>l) = 0) and therefore increase the density there. 
But Fig. [ll] shows a decrease in p up to cos(#l) w 0.2 and 
an increase towards smaller angles. Hence the gains in L z 
are big enough in order to deplete the polar regions and 
to tilt the orbital planes of the stars after their ejection 
to the equatorial plane. Moving in its plane, 8 pos of a star 
changes in the range 90° - p i anc < 6> pos < 90° + 6> p i anc , 
consequently resulting in a maximum of the density at 
= 90° (Fig. " 



10). So the EP changes its symmetry 



from spherical to cylindrical after the interaction with the 
binary. 

The stars which gain the most energy also gain the 
most angular momentum (see eq. (B.4)). In Fig. 12 we 
plot the energy versus cos(#l) = L z /L after the ejec- 
tion of stars. The figure shows a strong correlation be- 
tween the energy and L z /L, confirming the prediction of 
Eq. ( |B.4| ) : The orbital planes of the stars with the highest 
(kinetic) energies are more concentrated towards the equa- 
torial plane of the binary. A similar correlation is found 

when the eccentricity e + 2EL 2 is plotted versus the 

orientation of the orbits, showing that the most eccentric 
trajectories (e w 5.5) are most aligned with the plane in 
which the black holes rotate. 

The strong changes in L and E of the ejected stars 
compared to the almost constant values for the BP show 
that the BHB is mainly influenced by the ejected stars. 
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E as function of orbit orientation at t c 
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in 5pc distance, 308km/s. 6 m ax is the maximum impact 
parameter and assumed to be 50 pc since here the index of 
the density powerlaw changes from 2 to 4 and the density 
drops quickly with increasing radius (Fig. [IJ). For a mean 
density of 10 3 M pc -3 around the maximum value of 3.4 x 
10 3 M pc~ 3 at ~ 3pc and a = Ut yp /v3 we then obtain 



t 



relax 



5 x 10 1 



yr 



3 /l0 3 M Q pc- 3 
178km/s/ V P 



This time scale can now be compared to the time if r i c 
which a star needs to spiral from the outer edge to the 
inner edge of the torus due to dynamical friction. In fact, 
the diffusion coefficient is identical to the deceleration of 
a test star owing to dynamical friction by the surrounding 
field stars (Binney & Tremaine 1987), 



-0.2 0.2 0.4 0.6 0.8 1 

Lz I 

L lw 

Fig. 12. The energy at £ e j of the EP is plotted versus the 
orientation of the angular momentum. The more energetic 
the stars are, the more they are concentrated to the equa- 
torial plane. 

Provided that initially there »*-p enough stars the EP. can 



dv 

dt 



e r f(x) _ ^L- Xi 

\/7T 



Again the velocity is assumed to be distributed according 
to a Maxwellian and X = vj\[2a. The total number den- 
sity of the stars, with each having the mass m, is denoted 
by n, and erf is the error-function. The decelerating force 
is acting tangentially and thus causes the test star, also 
of mass m, to lose angular momentum per unit mass L 
at a rate dL/dt = rdv/dt. Following Binney & Tremaine 



gY t.rar J sufficient angular momentum, to »llnw the black ( lM7 D the star continues to orbit at the circular speed v c 



holes to merge, as shown in Sect. 3.2 



3.4. Does the torus survive the past-merger time? 

Finally we want to give in this section a crude estimate 
of the stability of the remaining torus-like distribution of 
the bound stars, once the BHs have merged. In a first 
approximation we can use the formulae from Sect. 2.1 for 



as it spirals to the center and its angular momentum per 
unit mass at radius r is at all times L = rv c . Equating the 
time derivative with the loss rate of L due to dynamical 
friction given above yields =r = ^%dt. Integrating this 
equation from a radius 10 pc, where the star starts to spi- 
ral inwards, to the inner edge at 2pc and solving for the 
time yields 



collisionless systems. With the following values, iVbn = 
AW(1 + N cj /N hn ) = 6.4 x 10 7 , R = 50 pc, 6 min = 2pc 
and the velocity set to the keplerian in a distance of 5 pc, 
v = 308km/s, we obtain for the ralaxation time 

iroiax = 4x 10 10 yr. 

In order to get a more precise time-scale we can make 
use of the diffusion coefficients D, which are a measure 
for the rate at which stars diffuse through phase spac e 
as a result of encounters, sec Binney fc Tremaine (1987 ). 
In the local approximation all encounters are assumed to 
be local, i.e. the impact parameter b is assumed to be 
much less than the system size R. Employing now the 
Fokker-Planck approximation and assuming a spherical 



kfric 
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The number density n(r) is approximated by the function 



n(f) 



,?.2 



1 



+ (W p i+(*r 



symmetry with a Maxwellian velocity distribution Binney 



fe Tremaine (1987) obtain for the relaxation time, defined 
as v 2 /D(Avt), 



t 



relax 



= 0.34- 



G 2 mp In A ' 

where A = b max v^ yp /2Gm. Here v typ is the typical velocity 
of stars of the system and set to the keplerian velocity 



which actually is a powerlaw with index —2 in the range 
f a < r < f j and index —4 for f > f&, while it vanishes for 
f < r a . 

We chose for the parameters r a — 2.75, f b — 50 and 
k = 5.3 x 10 4 pc -3 so that the volume integration from 
r = to infinity gives the right number of bound stars. 
Substituted in the expression for tf r j c and using for A, 
v c = ftyp and a the same values as above we obtain for 
the time a star needs to spiral from 10 pc distance to the 
inner edge of the torus at 2 pc 

i fric = 6.2 x 10 12 yr, 

which is in good agreement with i re i ax above. Thus we can 
conclude that the torus-like structure is stable. 



16 



C. Zier and P. L. Biermann: Binary black holes and tori in AGN 



Density /5(f) at the end t n 




le-10 



le-12 



0.1 



10 



100 



1000 



Fig. 13. At the end of the simulation the more stars are 
ejected the bigger M 2 is. While for q = 1 there is almost 
no cusp left in the center, the region close to the orbit of 
M 2 at f = 1 is not much depleted for q = 100. In the 
range 10 < f < 50 all density distributions follow the 
initial powerlaw with index —2 while at bigger distances 
p oc f~ 4 independent of q. The initial distribution is for 
all mass-ratios the same and represented by the solid line. 



4. Dependency on the mass-ratio 

Finally we will discuss the influence of the mass-ratio of 
the black holes on the ejected and bound star populations. 
In altering the mass-ratio q = Mi /Mi we always keep the 
mass of the primary black hole fixed to 1O 8 M . For the 
mass of the secondary we chose Mi = 10 8 and 10 6 solar 
masses corresponding to the ratios q — 1, 100 respectively. 
Hence also the total mass of the binary M m = Mi + Mi 
is changed. 

As Mi increases and q decreases the sphere of influence 
of the secondary BH expands and consequently the inner 
regions become more unstable. Thus a bigger fraction of 
stars is ejected and the density cusp bound to the primary 
BH only becomes less dense, being made up by 0.4% and 
2.8% for q = 1 and 100 respectively. Between r = 10r pc 
and 50 r pc the density is not much affected by variations 
in q and always follows a powerlaw with index —2 (see 
Fig. [ITJ). In the heated region (f > 50r pc ) too the density 
does not depend on the mass ratio and the relation p oc 
f -4 holds for all q. Due to the enhanced instability of the 
central region for a more massive secondary black hole 
the bound stars have to be more tangentially anisotropic 
in order to diminish the eccentricity of their orbits and not 
to come too close to the binary. This causes the inner edge 
of the torus at f ~ 2 r pc to be much sharper defined and 
thus the maximum of the density at f « 3 r pc to be much 
more pronounced for q = 1 , while it is almost smeared out 



for q = 100 (Fig. [ITJ). This is clearly visible in Fig. jhll 
where the density distribution is displayed at t = 10^^ 
and t = 10 7 73 yr and corresponds within a factor of a few 
to the lower limits of the merging times t = 10 6 - 89 and 
t = 10 7,77 yr of the BHs for q = 1 and 100 respectively. In 
fact, the bottom row for q = 100 shows no torus, but the 
density is diminished around the orbit of Mi and weakly 
indicates the formation of a shell-like distribution. This 
is confirmed if we let the simulation continue till £ max , 
about 10 times longer. Then the density distribution is 
comparable to that of q = 10 after 10 657 yr (2nd row in 
Fig-!). 

Thus the central stellar cluster in the potential of two 
merging BHs does not evolve into a torus-like structure 
for a sufficiently large mass-ratio q and rather stalls in the 
phase when it assumes a shell-like structure. 

On the other hand we can see that in the case of q = 1 
the final torus (t = 3.8 x 10 7 yr) is much more pronounced 
than for q = 10, compare Fig. || and [j~4| This means that 
for sufficiently small mass-ratios q a stellar torus forms 
in the late stages of a merger, while for young mergers 
and mergers with a big mass-ratio a shell-like density 
distribution is maintained. For q = 1 in Fig. [TJ| the in- 
ner edge of the torus is more sharply defined, and there 
is actually no cusp left in the center. Consequently the 
bound stars must have bigger angular momenta as Mi in- 
creases. If we compare the initial angular momentum of 
the ejected stars for the mass-ratios q = 1 and 100 with 
the ratio q — 10, it turns out that it increases with the 
mass of the secondary BH (decreasing with q) . The rela- 
tion (L(q = 1)) > (L(q = 10)) > (L(q = 100)) also holds 
at the end of the simulation as well as for the EP. Plotting 
the initial angular momentum versus the radius for each 
star for q = 1 and 100, just in the same way as for q — 10 
in Fig. ^, gives a very similar separation between the EP 
and BP. But with increasing q the transition value of L, 
dividing both populations, decreases from a little more 
than 2 for q = 1 to a value below 2 for q — 100 with 
always the same initial conditions. As the value of transi- 
tion of L in Fig. ^ is decreasing with Mi, also the mean 
values of the angular momentum of both populations are 
decreasing. Stars which have been ejected for small mass- 
ratios q "defect" to the BP as q increases and the inner 
region becomes more stable. Therefore the deviation AL 
of the BP increases while it is decreasing for the ejected 
population (see Table [l]) . 

The smaller the mass-ratio is, the more angular mo- 
mentum is extracted from the binary by a single ejected 
star on average ((L z (t ma , x )) — (L z (t = 0)), see Table |Tj). But 
this is not enough to compensate for the bigger amount 
of angular momentum of the binary due to the more mas- 
sive secondary BH, L m = L x + L 2 = 10 8 M 8 /(1 + q). Thus 
the cluster has to be more massive so that the binary 
can transfer its angular momentum to a larger number of 
ejected stars in order to merge. To compute the required 
number of stars and the time needed to enable the merging 
of th e bla ck holes we proceed in the same way as before in 
Sect. 3.2 for q — 10. With decreasing mass-ratio the max- 



C. Zier and P. L. Biermann: Binary black holes and tori in AGN 



100 



17 
100 




x=0 plane z=0 plane 

Fig. 14. Same as Fig. |^ for the final density distributions of the mass-ratios 5 = 1 (top row) and 100 (bottom row). 
For q = 1 there is no central cusp left and a pronounced torus with sharp defined edges emerges after 3.8 x 10 7 yr. 
There is no torus seen in case of q — 100 after 5.4 x 10 7 yr, only close to the orbit of the secondary BH the sensity is 
diminished. 



imum of the angular momentum lossrate is more peaked 
and has a value ~ 3.7 times higher for q — 1 compared to 
q = 10 (see Fig. [l5|). The maximum is reached the sooner, 
the smaller q is (t = 94, 142 corresponding to 1.34 x 10 5 
and 2.02 x 10 5 yr for q = 1, 10 respectively). Also the 
powerlaw with index —3/2 is approximated earlier by the 
loss-rate for small mass-ratios and therefore the system is 
relaxing faster. 

The number of stars and the times required for the 
merging are listed in Table [l]. The values we computed 



for the case q = 100 are not as accurate as for the other 
mass-ratios, since the data of the angular momentum loss- 
rate do not allow for a fit as good as for q = 1 and 10. 
However, the values show that for a binary with a fixed 
mass for the primary BH more stars in the cluster are re- 
quired for smaller mass-ratios in order to allow the black 
holes to coalesce. On the other hand the binary with the 
smaller mass-ratio is merging faster, provided that the 
cluster is sufficiently massive. For instance if q = 1, the 
merger is 2.4 times faster than the one where q = 10, while 
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Fig. 15. The smaller the mass-ratio q is, the more peaked 
is the maximum and the earlier the maximum is reached. 
Thus the system is relaxing faster. For all q the loss-rates 
finally approximate a powerlaw with index —3/2. For q = 
100 the curve is not as smooth due to the smaller number 
statistics. 

a factor of 3.5 times more stars in total are required. For 
a binary with smaller q the distance at r , where gravita- 
tional radiation starts to dominate the merging process, 
is bigger and thus reached after a shorter time due to 
the larger mass of M 2 (a tr « 1/120, 1/200, 1/265 pc and 
*t r » 6.18 x 10 6 , 1.40 x 10 7 , 4.71 x 10 7 yr for q = 1, 10, 100 
respectively, see Table |l|). The ratio of the time-ranges 
during which the merger is dominated by ejection of stars 
or the emission gravitational radiation, tt T /t mg , is indepen- 
dent of the mass-ratio and always equal to 4, see Eq. 

Assuming the average stellar mass to be that of the 
sun, these numbers show that the ejected mass M e j is of 
order of the total mass of the binary M m . While M j is 
decreasing more strongly with increasing q than the bi- 
nary mass, the ejected mass always exceeds that of the 
secondary BH. The ratio of both masses, M^jMi is in- 
creasing with q. Therefore a primary BH with mass M\ 
ejects more mass if it merges with TV black holes of mass 
Mi/TV than merging with another BH of mass Mi. This 



is the same result which has been obtained by Quinlan 
(1996|)]m scattering experiments. 

5. Conclusions 

One basic problem in our understanding of the central ac- 
tivity of galaxies is addressed in this paper: Do the black 
holes in the center of two merging galaxies also merge? 
Based on the two assumptions that galaxies contain su- 
permassive black holes in their centers and that galax- 
ies merge we simulated a stellar cluster in the potential 



of a BHB. When the distance between the black holes 
has shrunk to ~ r cusp and ejection of stars dominates the 
merging, we start our calculations. Since the simulated 
time is much smaller than the relaxation time of the clus- 
ter and the potential is dominated by the mass of the 
two black holes we were dealing with the restricted three 
body problem. The results we obtained by solving these 
equations for the stars of the cluster, initially distributed 
according to the singular isothermal sphere, just give an 
order of magnitude estimate since in the calculations we 
neglected further shrinking of the separation of the BHs 
due to ejection of stars. 

We find that the angular momentum divides the stellar 
cluster into a bound and ejected population. Stars staying 
bound in the potential of the BHs are distributed in a 
torus-like shape. Its relaxation time exceeds 10 10_12 yr, 
so that this structure can be regarded as stable. It can 
explain lots of the features which are postulated by the 
unification model for a dusty torus in AGN and will be 
the topic of paper 2. 

Stars which belong to the ejected population have low 
angular momentum and are moving on highly eccentric 
orbits, being radially anisotropic in velocity. The distri- 
bution of their peri- and apocenter always allows them 
to come very close to the binary. Here they gain suffi- 
cient energy in violent interactions with the black holes to 
become ejected. The ejected stars gain angular momen- 
tum in their z-component at the expense of the binary, 
which can be understood in terms of the conservation of 
the Jacobian Integral (dE oc dL z ). Thus the binary con- 
tinues to shrink. Integrating its loss-rate dL z /dt allows to 
calculate the number of ejected stars and the time required 
for the black holes to coalesce. 

Starting with lpc distance between the BHs, we find 
that a stellar cluster with a mass of order of the binary is 
required to allow the BHs to coalesce. The merging pro- 
ceeds on timescales of order of 10 7 yr and will chang e 
the spin of the central BH ( Wilson fc Colbert 1995 ). 
Consequently also the jet, emanating from the center, 
will jump and flow along the new spin direction after the 
merger is completed. This might be an explanation for the 
observed X-shaped radio galaxies, where a spectral aging 
time of the the secondary lobes suggests that a merger has 
happened ~ 6 x 10 7 yr ago (priv. comm. H. Rottmann). 
The agreement with our merger timescale is in favour of 
this idea. 

Most of the time (factor 4) this merging process is 
dominated by ejection of stars before emission of gravi- 
tational radiation takes over. Thus a stellar cluster with 
a total mass comparable to that of the binary is needed 
to allow the black holes to merge due to the ejection of a 
fraction of ~ 0.4 of all stars. After the ejection the stars 
are strongly concentrated to the equatorial plane of the 
binary, the more the higher their energy is, and the initial 
spherical geometry of the EP has been transformed into a 
cylindrical geometry. 

As the mass of the secondary black hole is increased, 
the central region becomes more unstable and a bigger 
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ejected stars bound stars 
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t 




(L) 


AL 


(L.) 


(L) 


AL 


iVbn 


•Ntotal 


ttt [yr] 


Otr [pc] 


(i) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


(11) 


(12) 


l 




^max 


0.086 
0.422 


1.165 
1.334 


0.532 
0.481 


-0.064 
-0.063 


3.773 
3.780 


1.432 
1.429 


0.667 


1Q 8.58 


1Q 6.79 


1/117 


10 




^max 


0.063 
0.284 


1.149 
1.264 


0.507 
0.497 


-0.045 
-0.043 


3.701 
3.711 


1.471 
1.467 


0.614 


10 7.97 


10 7.15 


1/197 


100 





0.086 
0.239 


1.005 
1.072 


0.411 
0.435 


-0.039 
-0.037 


3.411 
3.418 


1.562 
1.561 


0.395 


10 7.35 


1Q 7.67 


1/265 



Table 1. The lines in this table separate the results for the different mass-ratios q = M1/M2, as indicated in the first 
column. The second column subdivides each line into two lines for the next 6 columns, with the upper and lower line 
referring to numbers obtained in the beginning and the end of the simulation respectively. The 3rd and 6th column 
show the mean angular momentum of the z-component, the 4th and 7th column the mean total angular momentum 
and the 5th and 8th column the standard deviation of the total angular momentum for the ejected and bound stars 
respectively. The next column contains the number ratio N^/N\, n of ejected to bound stars and the 10th column the 
total number of stars which initially is required so that the ejected fraction extracts the complete angular momentum 
from the binary. i tr is the time needed by the black holes to shrink from an initial distance of 1 pc due to ejection of 
stars to atr, where emission of gravitational radiation starts to dominate further merging. 



fraction of stars is ejected. In order to extract sufficient an- 
gular momentum from the binary, more stars are required, 
so that the clusters mass still is of the order of magni- 
tude of the binary's mass. The system is relaxing faster 
and the distance where gravitational radiation dominates 
the merging process is reached earlier, so that the binary 
merges on smaller time scales (t w 10 6 ' 8 , 10 7 ' 2 and 10 7 ' 8 yr 
for q = 1, 10 and 100 respectively). The ratio M^jM^ is 
increasing with q and therefore a primary black hole with 
mass Mi ejects more mass in N minor mergers with sec- 
ondary black holes of mass M2 — M\/N than in one major 
merger with mass M2 — M\. 

In major mergers the mass distribution is likely to 
be completely rearranged and the resulting galaxy will 
probably have a rather elliptical than spiral shape. In 



their o verview of the s tructure of nearby AGN ( Malkan 
Gorjia: i, fc Tarn 1996 ) find that Seyfert 1 nuclei reside 
in more early- type host galaxies than Seyfert 2. This is in 
agreement with our finding, that the more massive the sec- 
ondary BH is, the more stars from the cluster are ejected. 
Thus the remaining torus-like distribution becomes shal- 
lower and there is a higher probability to see the Type 
1 nucleus directly. Combinig this with the conclusions by 



(Wilson & Colbert 1995), that radio loud objects are pos- 
sibly the product of major mergers, it is expected that 
radio loud galaxies are more likely to host Type 1 nuclei. 

The key result of this paper is: Provided that a black 
hole binary is surrounded by a cluster of solar mass stars 
with M c i uster w M m , i.e. the total mass of the cluster is 



comparable to that of the binary, the black holes merge on 
time scales of about 10 7 yr. In paper 2 we will demonstrate 
that the stellar torus surrounding the merged binary also 
meets the other requirements of the unification scheme. 
We will demonstrate that the stellar winds are pushed by 
radiation from the central source into long tails, and that 
the ensemble of many stellar wind-tails can constitute the 
observed torus. 
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Appendix A: Equations of the restricted three 
body system 

In the restricted three body system two massive bodies 
(BHs) are circulating around each other. A third parti- 
cle (a star), whose mass is negligible (m <C Mi,Mz), is 
moving in the time dependent gravitational potential of 
the massive bodies, which reads in the dimensionless form 

(eq. (§): $ = (^jr^ + . The star is located 

at f while f 1 and f 2 denote the positions of Mi and M2 re- 
spectively. The origin is identified with the center of mass 
of the two massive bodies and their axis of rotation is 
aligned with the z-axis and at t = with both black holes 
situated on the x-axis. The equations of motion for the 
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star are 



f = -V$ . 



(A.l) 



We always assumed Mi > M 2 (and therefore q = 
M\jMi > 1). The separation of the black holes, a = 
|f i — fa | = (1 + <z)|?i|, is constant, i.e. the black holes are 
moving on circular orbits with a constant radius around 
each other. The numerical integration of the equations of 
motion is processed faster in the comoving frame where 
the BHs are stationary on the i-axis in the distances x\ 
and X2 from the origin. The transformation into the co- 
moving frame introduces the centrifugal and the Coriolis 
force in (A.l): 



1 



(i x (ii x f ) 
f — f i 



26} X f 

f + qh 



fil 3 



i + q 

— ui x (cD x f ) — 2ui x f , 



gf i| 3 



(A.2) 



where u) in dimensionless units is actually the unit- vector 
in z direction, e z . Written in components with the vectors 
f = (x,y, z) and f = (v x , v y , v z ) this equation is trans- 
formed into a system of six coupled ordinary differential 
equations of first order: 



x = v x 

y = *v 

z = v z 

v x = - 







X 


-x, 


X 


f qxi 


l+« 


(■ 


f 




\r- 


-gfi| 3 




(■ 




^ 4 




V 


l+« 




r 


-fi| 3 


\r- 


- gfi| 3 




(« 




z 




z 


l+« 




r 


-fil 3 ^ 


\r- 


- gfi| 3 



+ 2v y + x (A. 3) 



2v x + y 



After being supplied with a set of initial conditions these 
equations have been integrated numerically. 

Appendix B: The Jacobian Integral 

For the restricted three body problem the energy is not 
conserved because the potential is time-dependent and 
consequently no integral of motion. The same holds for 
the angular momentum, which is not conserved since the 
potential is not invariant to rotations. But an integral of 
motion can still be found if we multiply the equation of 
motion for a star flA.2|) with f , 



(B.l) 



ff = — fV$ — (<D x x f) 

d d \ - Id,. _. 2 

— H — = ) $ H =(u> x tY 

dt dtj 2dV ' 

— (u> x r){ijj x f) . 

Because uj is constant and <fr does not depend explicitly 
on time we can rewrite this expression in a form showing 
that the quantity 

jE$ + if 2 -^xf) 2 , (B.2) 



the Jacobian Integral, is conserved (^ = 0). With the 
expressions for the energy and the specific angular mo- 
mentum expressed in the comoving frame 



1 f-2 

E = — [r + 2f (tD x f) + (u x f 



$(f) , 

L = rxf + fx(u>xf) 

and using the relation 

o>L = f (u; xf) + (iix f) 2 

we finally can write the Jacobian Integral in the simple 
form 



J = E — tDL = const. 



(B.3) 



Thus, in a rotating non-axisymmetric potential, neither 
E nor L is conserved, but the combination J = E — cDL. 



Applying infinitesimal variations to Eq. (BJ3), with £j = 
e z , leads to the expression 



dE = dL, 



(B.4) 



Consequently, stars which gain energy also gain angular 
momentum in the z-component. 



Appendix C: The isothermal sphere 

In a spherically symmetric potential there exist four iso- 
lating integrals, the energy E and the three components 
of the angular momentum L. According to the Jeans 
Theorem applied to spherical systems any non-negative 
function of these integrals can serve as the distribution 
function of a spherical stellar system. The extension of the 



Strong Jeans Theorem to spherical systems by Lynden- 



Bell (1962 ) allows to conclude that the distribution func- 
tion of any steady state spherical system can be expressed 
as function a f(E, L). If the system is spherically symmet- 
ric in all its properties, / is independent on the direction of 
L so that / = f(E, L). If now the potential $ is provided 
by the stellar system itself and we understand / as the 
mass distributionfunction (i.e. p = J fd 3 v) the Poisson 
equation 



A$ = 4nGp = 4nG / /d 3 v , 



(C.l) 



is the fundamental equation which governs spherical equi- 
librium stellar systems. In spherical symmetry this equa- 
tion is 



r 2 dr \ dr J 



4nG J f(^v 2 + <f>, |r x v|)d 3 v. (C.2) 



Tremaine (1987| ) 



Following the notation of Binncy 
we define the relative potential and the relative energy by 
^ = -$ + $o and £ = -E + $ = * - \v 2 . The constant 
<&o is choosen in a way that / > for £ > and / = 
for £ < 0. As <I>, the relative potential ^> also satisfies 
the Poisson equation = —AnGp with the boundary 
condition ^ — > $ f° r l r l ~~ * 00 ■ A distribution function 
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depending on the relative energy £ only further simplifies 
the spherical model. In such systems with / = /(\& — \v 2 ) 
the velocity dispersion tensor is isotropic everywhere be- 
cause the velocity dispersions are the same in all spherical 
components r, <p and 0. Using this distribution function 
and if we choose $o m a wa Y that f(£) = for £ < 0, 
Eq. (C.2) becomes 



r 



LA 



dr 



• 2 — = -16tt 2 G 



V27T 



'j 2 dv 



= -16tt 2 G / f(£)y/2(V-£)d£. (C.3) 
Jo 

Either we may understand this equation as a non-linear 
equation for \l/(r) with given distribution function /(£), 
or as a linear equation for f(£) with given ^>(r). With the 
choice for the distribution function for a stellar dynamical 
system 



/(*) = 



Pi 



(27TCT 2 ); 



Pi 



(2tt(7 2 ) ; 



■ exp 



±w 2 

2 v 



Poissons equation can be rewritten as 
d 

dr \ dr J 



_ (r 2 ^U 47 rGp. 



(C.4) 



(C.5) 



Integrating now the distribution function (C.4) over all 



velocities yields the density p = p 1 e*/ <T . With the help of 



this expression Eq. (C.5) becomes 



dr 



, dlnp 
dr 



AttG 



r 2 p . 



(C.6) 



Comparing this with the equation of hydrostatic equilib- 
rium for an inviscid fluid or gas 



d 
dr 



,dliip 
dr 



Gm 



4irr p 



(C.7) 



shows that they become identical if we choose the pa- 
rameter a to be er = hsT/m. This correspondence be- 
tween the gaseous and stellar-dynamical isothermal sphere 
originates in the velocity distribution of both systems. 
Integrating Eq. (CA) over the volume the distribution of 
velocities results in the Maxwell-Boltzmann distribution: 



F(v) = Ne 



2 /o 



(C.8) 



In order to find now a solution of Eq. (C.6) we simply 
assume a powe rlaw for the density (p = Cr~ b ) and substi- 
tute it in Eq. ( |06|) and find b = ^-Cr 2 - h . Thus the two 
parameters must have the values 6 = 2 and C = a 2 /2irG 
to fulfil this equation and we obtain for the density 



p(r) 



2irGr 2 



(C.9) 



A system described by this function for the density is 
known as the singular isothermal sphere, which is only 
true for big enough radii because of its singularity at the 
center. 
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